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Detailed  investigations  are  performed  to  analyze  and  mitigate  the  stability  issues 
encountered  in  developing  a  reduced  order  model  (ROM)  for  combustion  response  to 
specified  excitations  using  the  Euler  equations.  The  ROM  is  obtained  by  employing 
Galerkin’s  method  to  reduce  the  high-order  PDEs  to  a  lower-order  ODE  system  by  means  of 
POD  eigen-bases.  Possible  solutions  of  the  ROM  stability  issues  by  changing  and/or  by 
scaling  the  equation  variables  are  discussed  following  suggestions  from  previous  Euler 
equations  studies.  However,  our  evaluations  using  the  linearized  Euler  equations  indicate 
that  spurious  unstable  modes  are  still  encountered  in  the  resulting  ROMs.  Different  mean 
flow  and  boundary  conditions  are  implemented  to  further  evaluate  the  ROMs,  which 
indicate  that  the  presence  of  upstream  propagating  characteristic  waves  play  an  important 
role  in  affecting  ROM  stability.  Increasing  the  added  artificial  dissipation  terms  is  proposed 
and  shown  to  be  an  effective  method  that  insures  that  the  ROMs  are  both  numerically  stable 
and  capable  of  accurately  reproducing  the  CFD  solutions. 


I.  Introduction 

Combustion  instability  is  a  complex  phenomenon  that  results  from  the  coupling  between  the  modes  of  heat  release 
and  acoustics.  In  practical  combustor  devices  the  complexity  is  greatly  amplified  by  turbulent,  compressible 
flow,  very  high  rates  of  heat  release,  and  complicated  geometries  and  acoustic  boundary  conditions.  Modern 
computational  capability  offers  the  potential  for  moving  beyond  the  empirically-based  design  analysis  of  the  past, 
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but  high-fidelity  simulations  of  full  scale  dynamics  for  engineering  analysis  are  still  out  of  reach.  However,  high 
fidelity  simulations  of  smaller  scale  domains  can  be  used  to  obtain  reduced  order  models  of  the  combustion  response 
that  can  accurately  describe  the  linear/nonlinear  coupling  between  acoustics  and  combustion  and  can  subsequently 
be  used  to  analyze  full-scale  configurations. 

The  model  reduction  techniques  have  been  devised  to  develop  numerically  stable  and  robust  reduced  order 
models  (ROM)  [1-3]  and  they  have  been  applied  to  non-reacting  flow  problems  including  flow  control  [4-6]  and 
unsteady  aeroelasticity  [7,  8].  Recent  studies  have  extended  ROMs  to  combustion  problems  [9,  10].  A  preliminary 
exploration  of  the  POD/Galerkin  technique  using  a  model  reaction-advection  scalar  equation  for  developing  valid 
reduced-order  models  was  performed  by  Huang  et  al.  [11],  which  assessed  the  capability  of  the  ROM  for  predicting 
responses  at  target  frequencies.  This  work  was  further  extended  to  establish  ROM  performance  for  the  Euler  system 
of  equations  [16]. 

A  major  focus  of  these  previous  studies  focus  on  identifying  and  resolving  robustness  and  stability  issues 
associated  with  the  ROM  development.  As  reported,  the  issues  can  come  from  the  inherent  lack  of  numerical 
stability  in  the  POD/Galerkin  method  itself  [12],  truncation  of  low-energy  dissipative  POD  modes  [1]  and 
simplifications  of  higher-order  equations  [13].  The  balanced  POD  technique  has  been  proposed  to  build  numerically 
stable  ROM  for  linear  systems  [2,  14].  Bergmann  et  al.  proposed  to  add  residuals  of  the  Navier-Stokes  equations  to 
account  for  the  absence  of  low-energy  dissipative  POD  modes  [1].  Moreover,  Lucia  et  al.  [15]  demonstrated  the 
effectiveness  of  constructing  stable  ROM  by  including  additional  artificial  dissipation  terms.  In  addition  to  including 
artificial  terms  in  building  ROM,  researchers  have  also  tried  to  resolve  the  issues  based  on  the  numerical  properties 
of  the  system  equations.  Barone  et  al  [17,  18]  proposed  to  stabilize  the  reduced  system  by  symmetrizing  the  higher- 
order  PDE  with  a  preconditioning  matrix.  Rowley  et  al.  also  pointed  out  that  defining  a  proper  inner  product  can  be 
important  when  dealing  with  model  reduction  of  the  Navier-Stokes  equations  [3].  For  aeroelastic  applications, 
Amsallen  and  Farhat  have  shown  the  advantages  of  using  the  descriptor  form  over  non-descriptor  form  of  the 
governing  equations  [19]. 

It  is  important  to  point  out  that  the  stability  issues  in  the  above  studies  were  generally  associated  with  ROMs  that 
were  built  using  a  partial  or  incomplete  set  of  POD  modes.  Indeed,  our  previous  scalar  equation  studies  have  shown 
that  spurious  unstable  modes  can  be  generated  if  an  incomplete  set  of  POD  modes  (100%  of  the  energy  content)  is 
used  to  build  the  ROM  [11].  Conversely,  it  was  also  shown  that,  for  the  scalar  equation,  using  a  complete  POD  set 
always  insured  a  stable  ROM  [11].  However,  for  the  vector  system  of  Euler  equations,  Huang  et  al.  identified  that 
stability  issues  in  the  ROM  can  arise  even  when  a  complete  set  of  POD  modes  is  included  [16].  Based  on  the 
knowledge  above,  in  this  paper,  the  stability  issues  that  arise  in  the  ROM  development  of  the  vector  system  of  Euler 
equations  [16]  are  investigated  for  possible  causes  and  solutions  to  avoid  spurious  unstable  modes.  Again,  we  note 
that  the  complete  set  of  POD  modes  is  used  in  this  study  to  build  the  ROM  instead  of  a  partial  set  as  commonly  done 
in  reduced  order  development. 

The  remainder  of  the  paper  is  organized  as  follows.  In  Section  II,  we  present  the  Euler  equations  with  a  modeled 
combustion  source  term  and  the  forcing  function  and  boundary  conditions  used.  In  addition,  we  also  present  the 
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Galerkin  formulation  and  the  POD  techniques  for  deriving  reduced  order  models  for  the  linearized  version  of  the 
Euler  equations.  In  Section  III,  a  summary  is  given  with  regard  to  previous  efforts  to  build  stable  ROMs  by  changing 
the  solution  variables  used  for  the  POD  generation  [16].  In  Section  IV,  we  describe  the  test  problem  and  summarize 
the  mean  flow  conditions  used  for  ROM  studies.  In  Section  V,  we  demonstrate  the  existence  of  stability  issues  and 
the  corresponding  consequences  in  ROM  development.  We  then  use  different  mean  flow  and  boundary  conditions  to 
evaluate  their  effects  on  ROM  stability  and  draw  some  conclusions  regarding  when  the  stability  problems  occur. 
Following  this,  we  propose  and  test  the  inclusion  of  additional  artificial  dissipation  to  mitigate  the  stability  issues.  In 
the  final  section,  we  provide  concluding  remarks  and  suggest  future  directions  for  continued  research. 


II.  Formulation 


A.  Governing  equations 

The  governing  equations  are  the  quasi-one-dimensional  unsteady  Euler  equations  with  a  single-step  chemical 
reaction  and  a  specified  reaction  distribution. 
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Elere  x  and  t  are  the  space  and  time  variables,  p  is  the  density,  u  is  the  velocity,  e  is  the  total  energy,  p  is  the 
pressure,  Yox  is  the  oxidizer  mass  fraction  and  A  =  A(x)  is  the  cross-section  area  of  the  geometry.  The  effects  of  fuel 
addition  are  accounted  through  steady  source  term  H f,  where  C0f  =  C  l(/box  with  constant  Cf!0  representing  fuel-to- 
oxidizer  ratio  and  a  sinusoidal  spatial  distribution  is  used  to  model  the  oxidizer  reaction. 
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and  end  of  the  combustion  zone.  The  reaction  constant  kf  is  selected  to  ensure  that  the  oxidizer  is  consumed  within 
the  specified  combustion  zone.  The  unsteady  combustion  response  is  accounted  in  the  source  term,  H(;,  using  the  n-r 
■w  ~  t  \  p{xP~r)~p{x) 

model  [20],  q  =  q-n-aix ) - — - ,  which  relates  the  unsteady  heat  release  to  the  pressure  oscillations 

P\x) 


through  an  index,  n  and  a  time  lag  constant,  r.  Here  q  is  the  integral  mean  combustion  heat  release  and  a(x)  is  a 
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,  v  1 

scaling  function  following  a  normal  distribution,  a[x)  =  — ^=exp 

G'IIk 


2(7  2 


.  This  model  was  previously  used 


[21]  to  simulate  combustion  instability  in  a  longitudinal  rocket  combustor, 

For  simplicity,  the  linearized  version  of  Eq.  (1)  is  used  for  the  studies  in  ROM  development, 
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B.  Boundary  conditions 

Fluctuating  conditions  inside  the  computational  domain  are  obtained  by  specifying  periodic  oscillations  of  target 
quantities  ),  which  can  be  mass  flow  rate  ( nV(t) ),  back  pressure  (  p'back )  or  the  Riemann  variables,  about  their 
mean  value  ( <j)Q ), 

Nf 

f(0  =  $)^-Zsin(27r(/o+4/>)  (3), 

A  f  k= i 

where  f0  is  the  initial  frequency;  A f  represents  the  fundamental  frequency  increment  and  Ay  is  the  total  number  of 
frequencies  included  in  the  forcing  function.  The  fundamental  period  Tp  is  determined  by  A f  such  that  Tp  =  1  /A/.’ 

Note  that  if  A f=fo,  Eq.  (3)  represents  a  standard  Fourier  series  although  herein  we  generally  take  A f  <  f0  so  that  Eq. 
(3)  differs  from  a  Fourier  series. 

Four  different  boundary  conditions  (summarized  in  Table  1)  are  implemented  to  help  further  understand  and 
identify  the  causes  of  numerical  stability  issues  arising  from  the  previous  ROM  study  for  the  Euler  equations  [16]. 


B.C.  label 

Upstream 

Downstream 

Subsonic  flow 

1 

m ,  T°  and  Yox 

PbdiCk. 

2 

Riemann  invariant  B.C. 

Non-reflecting 

3 

m ,  T°  and  Yox 

Non-reflecting 

Supersonic  flow 

4 

p° ,  m ,  T°  and  Yox 

/ 

Table  1  Summary  of  boundary  conditions. 
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C.  Construction  of  POD  eigen-bases  for  vector  equations 

POD  eigen-bases  are  calculated  based  on  the  CFD  solutions,  Q'p  ( x,  t  j  obtained  from  Eq.  (2)  using  the  vector¬ 
valued  method. 
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where  a„  is  the  singular  value  of  the  nth  POD  mode  (scalar),  a„(1)  is  the  nth  POD  temporal  mode  (scalar),  and  the  nth 
eigen-mode,  On(x),  is  an  orthonormal  vector  function. 


1,  if  k  =  n 
0,  otherwise’ 


X  =  {0  <x<L) 


(5). 


It  should  be  noted  that,  unlike  in  the  scalar  equation  case,  to  obtain  reasonably  scaled  POD  eigen-bases  for  the 
vector  variables,  a  normalization  matrix  P(x)  must  be  used  before  calculating  the  POD  eigen-bases.  Likewise,  to 
reconstruct  the  CFD  solutions  using  the  POD  eigen-bases,  the  matrix  P(x)  again  needs  to  be  included, 


Q  'Pix’t)apl(x)llan(t) 


t,n  (*) 
<t>T,n  (*) 

K.n  (X) 


(6). 


Different  definitions  of  P(x)  are  given  later. 


D.  Model  reduction  of  Euler  equations 

The  application  of  the  POD-Galerkin  method  to  the  linearized  Euler  equations  (Eq.  (2))  is  briefly  introduced  here. 
Additional  details  can  be  found  in  Ref  [11].  Upon  obtaining  the  eigen-bases  as  in  Eq.  (4),  the  target  governing 
equation  is  projected  onto  the  kth  eigen-mode,  0/(x),  throughout  the  whole  computational  domain.  Before  the 
projection  the  governing  equation  needs  to  be  normalized  by  pre-multiplying  by  the  matrix  P(x), 
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(7), 


Substituting  the  POD  expansion,  Eq.  (6)  into  Eq.  (7)  and  using  a  numerical  quadrature  to  approximate  the  integrals, 
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Following  the  model  reduction  procedure  in  Ref  [11]  and  using  a  consistent  discretization  scheme  (a  2nd-order 
upwind  scheme  is  used  in  both  CFD  and  ROM)  to  approximate  the  gradient  term  in  Eq.  (8),  an  ODE  system  can  be 
obtained  with  the  contributions  from  boundary  conditions  appearing  as  a  source  term  on  the  right-hand-side, 


IMl)  h(/) 


(9), 


where  a (f)  =  [a,(0  -  ak(t)  ■■■  aNf  (t)J , 

h(r)  =  l^/j,  (r)  hk  (r)  •••  hN  (r)J  with  contributions  from  boundary  conditions,  hk  (?)  =  Fk  (^'(f)) , 

L  is  the  stiffness  matrix  which  describes  the  dynamics  of  the  reduced  ODE  system. 


III.  Background 

Our  previous  studies  with  the  scalar  equation  have  demonstrated  that  using  a  complete  set  of  POD  modes  is 
sufficient  to  insure  stability  in  the  ROM  development  [11].  However,  perhaps  surprisingly,  our  studies  with  the 
vector  system  of  Euler  equations  indicate  that  this  is  no  longer  the  case  [16].  In  other  words,  spurious  unstable 
modes  are  encountered  even  when  the  ROM  Is  developed  using  the  complete  set  of  POD  modes.  In  fact,  this 
problem  arises  even  for  the  simple  problem  of  acoustic  modes  in  a  straight  duct  without  combustion,  which  should 
be  stable  by  definition.  It  has  also  been  demonstrated  in  our  previous  study  that  different  scaling  strategies  of  the 
solution  variables  can  be  helpful  in  getting  rid  of  the  spurious  unstable  modes  [16].  Moreover,  as  mentioned  in  the 
introduction,  it  has  been  reported  by  other  researchers  that  an  appropriate  treatment  of  solution  variable  for  the  POD 
mode  calculation  can  be  useful  in  resolving  the  ROM  stability  issues.  Therefore,  the  effects  of  different  solution 
variable  treatments  on  the  ROM  stability  are  briefly  discussed  in  this  section. 

First,  the  scaling  of  temperature  fluctuations  (7”)  is  tuned  through  parameter,  a,  by  defining  the  matrix  P(x)  in 
Eq.  (4)  as, 

P{*)  =  diag  (1  /  PL c  .  1  /  «C< , «  /  C .  1  /  C,™*  )  (10), 

where  p'max  =  Max  (  //(x,z)} ,  V  0  <  x  <  L  and  t  so  that  the  variations  in  each  variable  can  be  taken  into  account  in 

relative  to  their  maximum  amplitude.  As  an  example,  CFD  solutions  are  generated  using  B.C.  1  in  Table  1  by 
perturbing  the  inlet  mass  flow.  Two  specific  forcing  functions  are  considered  with  3  and  5  frequencies  included.  As 
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in  the  previous  studies  [16],  the  eigen-values  of  stiffness  matrix  L  in  Eq.  (9)  are  investigated  to  assess  the  numerical 
stability  of  the  resulting  ROMs.  The  real  parts  (cr)  of  the  eigen-values  are  plotted  against  the  scaling  parameter,  a, 
for  the  two  forcing  functions  in  Fig.  1.  It  is  noted  that  a  stable  ROM  (cr<  0)  can  be  achieved  for  an  appropriate  value 
of  the  scaling  parameter.  However,  the  two  cases  shown  here  do  not  share  a  common  scaling  strategy  that  insures 
stability  in  both  ROMs,  which  means  that  the  choice  of  scaling  is  problem-dependent. 


3 -frequency  forcing 
(1000,  1250  and  1500Hz) 

Figure  1  The  scaling  effects  on  the  eigen-values  < 
variables  (the  stable  region,  a 


5-frequency  forcing 
(750,  1000,  1250,  1500  and  1750Hz) 

resulting  ROMs  using p',u’,  T’  and  Y’ox  as  solution 
0,  is  covered  by  the  grey  box). 


As  a  second  example,  different  solution  variables  (summarized  in  Table  2)  from  the  same  CFD  simulation  are 
selected  to  build  the  ROM.  The  CFD  solutions  are  generated  using  the  3-frequency  forcing.  The  set  1  in  Table  2 
uses  the  same  solution  variables  {p  \  u’,  T’  and  Y’ox)  as  in  Fig.  1  and  all  variables  are  scaled  to  the  same  fluctuating 
levels  (a=  1  in  Eq.  (10)).  Set  2  uses  the  same  scaling  strategy  as  set  1  but  instead  of  using  temperature  fluctuations 
(7”),  density  ( p ’)  is  used.  Set  3  follows  the  strategy  derived  by  Barone  et  al.  [17].  Their  method  symmetrizes  the 
governing  equations  by  pre-multiplying  a  matrix  S  and  choose  the  solution  variables  as 

Qs =(«'  (i ip)  p  ^:)T , 


"p  0  0  (P 

0  ypp 1  p  0 

0  p  0 

YP 

v  0  0  0  YP; 


(11), 


7 


where  /is  the  ratio  of  the  constant  volume  and  constant  pressure  heat  capacity.  And  then  the  POD  modes,  CD,,. ,  are 

Np 

calculated  based  upon  the  transformed  variables  Q'T  =  S1  ~  ^ak  (f)<D(.  (x) . 

k= 1 


Case  Name 

Variables 

Set  1 

p\u  \  T’  and  Y’ox 

Set  2 

p\u\p  and  Y’ox 

Set  3 

Symmetrizing  variables  [17] 

Table  2  Summary  of  solution  variables  used  for  ROM  construction. 


The  ROM  spectrum  comparisons  are  shown  in  Fig.  2  for  the  three  solution  variables.  The  real  (cr)  and  imaginary 
(271 f)  parts  of  the  stiffness  matrix  L  eigen-values  are  plotted  against  each  other.  Unstable  eigen-modes  (cr  >  0, 
highlighted  using  dashed  circle)  can  be  seen  using  solution  variables  set  1  and  3  while  a  stable  ROM  can  be  obtained 
using  set  2  variables.  On  the  other  hand,  when  different  forcing  functions  are  used,  we  find  that  ROM  stability 
cannot  be  guaranteed  even  for  set  2.  Moreover,  the  symmetrizing  variables  in  set  3,  which  was  reported  to  be 
capable  of  generating  stable  ROMs,  also  seem  to  have  stability  issues  here.  In  addition,  we  have  also  tried  the 
conservative  set  of  variables  (p\  ( pit )  ’,  ( ph°-p )  ’  and  (pYox)  ’)  following  Amsallen  and  Farhat’s  studies  [19]  but  again, 
we  find  that  ROM  stability  cannot  be  universally  guaranteed. 
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Figure  2  The  ROM  spectra  comparisons  for  different  selections  of  solution  variables  for  the  3-frequency 

forcing  case  (1000, 1250  and  1500Hz). 


In  summary,  the  ways  to  generate  POD  modes,  in  terms  of  either  scaling  of  vector  variables  or  the  selection  of 
solution  variables  can  be  important  in  affecting  the  stability  properties  of  resulting  ROMs.  However,  it  seems  to  be  a 
difficult  task  to  look  for  a  universally  valid  method  for  POD  mode  calculation  to  resolve  the  ROM  stability  issues. 
Therefore,  a  further  systematic  investigation  is  needed  to  provide  insights  regarding  possible  solutions  of  the 
stability  issues  in  reduced  order  model  development. 


IV.  Overview  of  Test  Problem 

A  constant-area  duct  is  set  up  for  the  ROM  development  as  shown  in  Fig.  3  with  different  boundary  conditions 
summarized  in  Table  1.  For  a  typical  computation,  the  steady  state  numerical  solution  of  Eq.  (1)  is  computed  first. 
The  upstream  boundary  condition  is  then  switched  to  a  temporally  varying  periodic  condition  as  defined  in  Eq.  (3) 
and  the  time -accurate  computation  is  continued  until  stationary  conditions  are  reached. 


x/L  =  0 


x/L  =  1 


Upstream 

B.C. 


Downstream 

B.C. 


Figure  3  Computational  setup  for  the  one-dimensional  Euler  equations. 


The  solutions  are  tabulated  and  stored  periodically  (i.e.,  at  several  time-levels  within  a  forcing  period)  thereby 
generating  a  rectangular  matrix  that  can  be  used  as  a  database  for  fitting  eigen-bases  by  means  of  the  POD  procedure 
in  Eq.  (4).  The  POD  eigen-bases  are  then  applied  to  the  governing  linear  or  nonlinear  PDE  to  derive  the  reduced- 
order  ODE  formulation  (or  ROM).  In  order  to  focus  our  attention  on  understanding  the  numerical  stability  issues 
identified  in  Ref  [16],  no  combustion  is  included  in  Eq.  (1)  and  (2)  and  only  gas  dynamics  (p’,u  ’  and  7”)  are  solved 
in  the  constant-area  duct.  The  five  mean  flow  conditions  simulated  correspond  to  a  range  of  Mach  numbers 
summarized  in  Table  3. 


Case  # 

Ma 

Pressure,  MPa 

T emperature,  K 

Subsonic 

A 

0.21 

1 

2000 

B 

0.46 

1 

1900 

C 

0.86 

1 

1670 

Supersonic 

D 

1.87 

4 

2000 

E 

2.80 

1.5 

2000 

Table  3  Summary  of  mean  flow  conditions. 


9 


V.  Results 


A.  Baseline  ROM  Characteristics 

Case  A  in  Table  3  is  selected  as  the  baseline  for  ROM  characteristics  studies  using  B.C.  1  in  Table  1.  It  has  been 
reported  in  previous  work  that  the  ROM  stability  can  vary  by  applying  different  forcing  functions  in  obtaining  CFD 
solutions  using  the  same  boundary  conditions  [16].  Therefore,  six  forcing  functions,  from  single  to  multiple 
frequencies,  are  used  for  ROM  generation  and  are  summarized  in  Table  4.  The  POD  eigen-basis  is  generated  using 
the  vector-valued  method  defined  in  Eq.  (4)  and  the  three  variables  ( u  ’,  p  ’,  and  T’)  are  normalized  by  the  maximum 
values  of  the  respective  fluctuating  quantity  throughout  the  whole  domain  (i.e.  p(x)  =  diag  (1  /  p'nmy ,  1  /  u'max ,  1  /  T'm  ) 

where  p'ma  =  Max  ,  V  0  <x<L  and  t  etc.).  Following  the  conclusion  from  scalar  Equation  studies  [11], 

all  the  meaningful  (or  non-zero  amplitude)  POD  eigen-bases  are  used  to  generate  the  ROM  corresponding  in  each  of 
the  six  cases. 


Case  Name 

fo,  Hz 

A/  Hz 

#  of  forcing  frequencies 

#  of  necessary  POD  modes 

Singe-frequency 

-- 

1 

2 

2-frequency 

250 

2 

4 

3-frequency 

1000 

250 

3 

6 

4-frequency 

250 

4 

8 

10-frequency 

250 

10 

18 

21 -frequency 

250 

21 

18 

Table  4  Summary  of  cases  used  in  studies  of  ROM  characteristics  for  Case  A  in  Table  3. 


The  eigen-value  spectra  of  the  ROM  stiffness  matrix  (L  in  Eq.  (9))  are  investigated  and  shown  in  Fig.  4  for  the 
six  cases  in  Table  4.  For  the  single-,  2-  and  3-frequency  forcing  cases,  all  eigen-values  are  observed  to  have  negative 
real  parts  (cr),  which  indicates  that  the  corresponding  ROMs  will  be  stable.  Solving  the  ROM  equation  for  these 
three  cases  using  the  same  forcing  as  used  in  the  CFD  solution  gives  reconstructed  solutions  that  are  indeed  in 
excellent  agreement  with  the  CFD  results  as  shown  in  Fig.  5  (top).  The  ROMs  are  able  to  reproduce  the  original 
CFD  solutions  and  show  no  evidence  of  unphysical  or  unstable  behavior. 
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Figure  4  Eigen-value  spectrum  of  ROM  stiffness  matrix  for  Case  A  in  Table  3. 


Single-frequency 


2-frequency 


3 -frequency 


4-frequency 
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0.028  0.03 
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Figure  5  ROM  and  CFD  solutions  comparison  of  Case  A  at  x/L  =  0.5  for  cases  in  Table  4. 
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The  remaining  three  cases  with  multiple  frequencies  in  the  forcing  function  appear  to  be  more  problematic  for 
building  stable  ROMs.  Unstable  modes  can  be  observed  in  the  4-,  10-  and  21 -frequency  cases  (highlighted  by  the 
dashed  circle  in  Fig.  4)  even  though  a  complete  set  of  POD  modes  have  been  included  in  building  the  ROM.  These 
unstable  modes  leads  to  unphysical  growth  in  the  ROM  solution  as  shown  in  Fig.  5  (bottom).  With  higher  growth 
rate  of  the  unstable  eigenvalue,  the  ROM  prediction  is  seen  to  go  off  track  faster.  In  fact,  our  results  suggest  that  the 
more  diverse  the  frequencies  are  in  the  forcing  function,  the  more  difficulties  the  ROM  construction  can  encounter. 


B.  Potential  causes  for  stability  issues  in  ROM 

To  further  understand  and  identify  the  potential  causes  for  the  unstable  modes  in  the  ROMs  observed  in  Fig.  5, 
two  major  aspects  are  investigated  in  this  section: 

1.  So  far,  the  studies  have  focused  on  low  Mach  number  flows  (Ma  <  0.3);  so,  we  study  the  effects  of  flow 
conditions— specifically  the  Mach  number  variation— on  the  ROM  stability  characteristics.  Importantly,  we 
consider  the  difference  between  subsonic  flow  that  have  characteristic  waves  traveling  in  both  directions  and 
supersonic  flow  where  the  waves  are  all  traveling  downstream.  Importantly,  we  note  that  the  latter  case  more  closely 
resembles  the  scalar  wave  equation  case. 

2.  Based  on  the  knowledge  from  previous  studies  [11],  we  know  that  the  stability  issues  do  not  arise  in  the  scalar 
equation  studies  as  long  as  a  complete  set  of  POD  modes  are  included  in  the  ROM  construction.  We  therefore  seek 
to  isolate  the  characteristic  waves  contained  in  the  Euler  equations  through  control  of  the  boundary  conditions  as 
given  in  Table  1.  Specifically,  the  Riemann  and  non-reflective  boundary  conditions  serve  to  insure  that  the 
characteristic  waves  stay  decoupled,  which  should  mimic  the  scalar  situation  pretty  closely.  On  the  other  hand,  the 
use  of  reflective  boundary  conditions  introduces  wave  reflections  and  interactions  that  are  unique  to  the  vector 
system. 

The  4-frequency  forcing  function  in  Table  4  is  selected  as  the  baseline  perturbation  to  obtain  CFD  solutions.  The 
POD  eigen-basis  is  generated  based  on  the  vector-valued  method  and  the  maximum  values  of  each  fluctuating 
quantities  are  used  for  normalization  as  mentioned  in  section  A. 


a.  Effects  of  Mach  number 

The  ROMs  are  generated  based  upon  different  conditions  summarized  in  Table  3  from  low  to  high  Mach 
number.  The  ROM  eigen-value  spectrum  for  subsonic  flow  cases  (A,  B  and  C)  is  shown  in  Fig.  6  and  it  can  be  seen 
that  for  low  Mach  number  (<0.5),  unstable  modes  (highlighted  in  dashed  circle)  are  still  present  while  when  the 
Mach  number  gets  higher  (>0.8),  the  resulting  ROM  becomes  stable  (negative  real  parts,  a).  The  reconstructed 
ROM  solutions  using  the  4-frequency  forcing  function  in  Table  4  are  compared  against  the  CFD  solutions  for  Case 
B  and  C  in  Fig.  7.  As  expected  and  similar  to  Fig.  5,  unphysical  amplitude  growth  that  is  markedly  different  from 
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the  original  CFD  solution  is  seen  for  Case  B,  while  a  stable  ROM  is  obtained  for  Case  C,  which  correspondingly 
reproduces  the  original  CFD  solutions  excellently. 


Figure  6  Eigen-value  spectrum  of  ROM  stiffness  matrix  for  Case  A,  B  and  C  (subsonic  flow)  in  Table  3  using 

4-frequency  forcing  function. 


Case  B  Case  C 

Figure  7  ROM  and  CFD  solutions  comparison  at  xlL  =  0.5  for  Case  B  (left)  and  C  (right)  in  Table  3  using  4- 

frequency  forcing. 


As  it  goes  to  supersonic  flow  cases  (D  and  E),  all  the  resulting  ROMs  are  stable  based  on  the  eigen-value 
spectrum  shown  in  Fig.  8  with  all  the  real  parts  of  eigen-values  being  well  below  zero.  The  observations  from  Figs. 
6  and  8  provide  some  clues  as  to  the  possible  causes  of  the  unstable  eigen-modes  in  ROMs.  As  the  Mach  number  of 
flow  goes  higher,  the  reflecting  characteristic  waves  (or  upstream  traveling  waves)  from  downstream  boundary  are 
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weakened  and,  when  the  flow  becomes  supersonic,  there  is  no  wave  traveling  back  from  the  downstream  end.  This 
suggests  that  the  ROM  stability  issues  might  be  related  to  the  simultaneous  presence  of  downstream  and  upstream 
traveling  waves  in  the  CFD  solution— a  situation  that  only  arises  in  the  subsonic  Euler  system  (and  not  in  supersonic 
flow  or  in  the  scalar  equation  case).  In  turn,  this  motivates  the  study  in  the  following  section  that  studies  the  effects 
of  boundary  conditions  which  are  also  directly  related  to  the  wave  characteristics  of  the  CFD  solution. 
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Figure  8  Eigen-value  spectrum  of  ROM  stiffness  matrix  for  Case  D  and  E  (supersonic  flow)  in  Table  3  using 

4-frequency  forcing  function. 
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b.  Effects  of  boundary  conditions 

Riemann  invariant  boundary  conditions  are  implemented  at  both  the  upstream  and  downstream  boundaries  in 
order  to  eliminate  wave  reflections.  The  perturbations  defined  in  Eq.  (3)  are  applied  for  each  of  the  characteristics 
waves  individually  using  the  4-frequency  forcing  function  in  Table  4.  By  doing  so,  the  downstream  and  upstream 
propagating  waves  are  separated  so  that  the  waves  do  not  interact.  The  Riemann  variables  correspond  to  propagating 
speeds,  u,  u+c  and  u-c.  Three  cases  are  considered  here:  1.  only  u  and  u+c  characteristic  waves  are  perturbed  while 
the  downstream  propagating  wave  (u-c)  is  kept  unperturbed;  2.  only  u  and  u-c  characteristic  waves  are  perturbed 
while  the  downstream  propagating  wave  (u+c)  is  kept  unperturbed;  and  3.  all  characteristic  waves  are  perturbed. 
The  eigen-value  spectrum  of  the  resulting  ROMs  of  these  three  cases  are  shown  in  Fig.  9.  It  can  be  readily  seen  that 
all  ROMs  are  stable  with  negative  real  parts,  a,  which  confirms  that  as  long  as  the  waves  are  treated  independent  of 
each  other,  the  Euler  system  reduces  to  a  set  of  three  scalar  wave  equations  (at  least  in  the  linear  limit)  and  the 
resulting  ROM  problem  remains  stable. 
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Figure  9  Eigen-value  spectrum  of  ROM  stiffness  matrix  for  B.C.  2  in  Table  1  using  4-frequency  forcing 

function. 


In  the  second  instance,  the  mass  flow  rate  is  specified  for  upstream  boundary  (and  forced),  while  a  non-reflecting 
boundary  condition  is  used  downstream  to  eliminate  wave  reflections  into  the  upstream  propagating  wave  (m-c).  The 
eigen-value  spectrum  of  the  resulting  ROM  for  this  case  is  shown  in  Fig.  10,  which  also  does  not  show  any  unstable 
modes.  This  result  should  be  contrasted  with  Case  A  in  Fig.  6,  which  specified  constant  back  pressure  for  the 
downstream  boundary  thereby  allowing  a  reflected  upstream-running  characteristic  wave.  It  can  therefore  be 
inferred  that  the  boundary  reflection  of  the  upstream  propagating  wave  (m-c)  is  a  potential  cause  for  the  ROM 
stability  issues. 

Based  on  the  results  above,  it  can  be  concluded  that  the  ROM  stability  issue  arises  when  both  considerable 
upstream  and  downstream  propagating  characteristics  waves  are  present  and  mixed  at  the  boundaries.  Moreover, 
stable  ROMs  can  be  obtained  in  two  major  ways:  (1)  by  reducing  the  upstream  propagating  wave  by  either 
increasing  the  flow  Mach  number  or  (2)  by  applying  non-reflecting  boundary  conditions  to  avoid  wave  reflections. 
These  findings  clearly  confirm  why  the  stability  issues  occur  in  the  vector  Euler  system  case  and  not  in  the  scalar 
equation  studies.  Of  course,  further  study  is  needed  to  understand  precisely  why  the  ROM  solutions  are  so 
influenced  by  the  presence  of  boundary  wave  reflections  and  will  be  undertaken  as  part  of  future  work.  For  now,  we 
examine  the  possibility  of  mitigating  the  stability  problems  by  the  addition  of  artificial  dissipation  terms.  This 
technique  was  originally  proposed  for  ROM  development  using  an  incomplete  POD  eigen-bases,  but  we  use  it  here 
for  situations  wherein  a  complete  eigen-bases  is  used  in  the  ROM  construction. 
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Figure  10  Eigen-value  spectrum  of  ROM  stiffness  matrix  for  B.C.  3  in  Table  1  using  4-frequency  forcing 

function. 


C.  ROM  stabilization  through  additional  artificial  dissipation 

It  has  been  shown  in  previous  studies  that  a  simple  treatment  like  including  additional  artificial  dissipation  in 
building  the  ROM  can  useful  to  eliminate  the  non-physical  unstable  modes  [15].  Here  the  method  of  including  the 
artificial  dissipation  is  briefly  introduced  starting  with  the  discretized  form  of  the  linearized  model  equation,  Eq.  (2), 
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where  Ax,  is  the  size  of  the  zth  cell  and  areai+in  is  the  left/right  faces  area  of  the  ;th  cell.  The  artificial  dissipation  is 
added  at  the  cell  faces, 
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(13). 


It  should  be  noted  that  the  artificial  dissipation  already  exists  in  the  original  CFD  solutions  since  we  are  using  a  2nd- 
order  upwind  discretization  scheme  (with  ft  =  1),  based  upon  which  the  POD  eigen-basis  is  calculated.  To  increase 
the  stabilizing  influence  of  the  ROM,  different  ft  values  (>1)  are  used  during  the  Galerkin  projection  step.  By 
systematically  varying  the  definition  of  the  /?  parameter,  we  can  estimate  how  much  extra  artificial  dissipation  is 
needed  to  stabilize  the  ROMs. 
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Case  A  (4-frequency) 


Case  A  (10-frequency) 


Figure  11  Eigen-value  spectrum  of  ROM  stiffness  matrix  by  including  different  amounts  of  artificial 

dissipation 


The  four  cases  showing  unstable  modes  in  the  resulting  ROMs  are  selected  from  section  A  and  B  for  further 
investigation.  These  are  Case  A  (4-,  10-  and  21-frequency)  and  Case  B  (4-frequency).  The  corresponding  ROM 
eigen-value  spectra  are  shown  in  Fig.  1 1  with  different  amounts  of  artificial  dissipation  added.  It  can  be  readily  seen 
from  the  spectra  that  including  additional  dissipation  (increasing  the  /?  value  greater  than  unity)  is  able  to  stabilize 
the  ROMs  by  reducing  the  real  part  (er)  of  each  mode  while  it  does  not  materially  affect  the  imaginary  part  (/)  very 
much.  Moreover,  it  is  also  observed  to  have  a  more  significant  effect  on  the  unstable  modes  than  on  the  stable 
modes,  which  is  a  desired  property  since  we  do  not  want  to  alter  the  physically  accurate  eigenmodes.  For  cases  with 
less  frequency  diversity  in  the  forcing  function  (4-frequency),  the  real  part  of  the  unstable  mode  is  reduced  below 
zero  when  /?  is  increased  up  to  2.2  ~  2.4,  which  means  approximately  twice  the  artificial  dissipation  is  needed  to 
generate  a  stable  ROM.  For  cases  with  more  frequency  diversity  in  the  forcing  function  (10-  and  21 -frequency), 
approximately  3  times  more  artificial  dissipation  is  needed  (/?-  3.7)  to  stabilize  the  ROMs. 
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Figure  12  ROM  and  CFD  solutions  comparison  at  xlL  =  0.5  by  including  additional  artificial  dissipations. 


The  stabilized  ROMs  with  the  corresponding  0  values  are  tested  and  validated  against  the  CFD  solutions  in  the 
same  way  as  in  Figs.  5  and  7.  The  reconstructed  ROM  solutions  are  shown  in  Fig.  12.  All  the  stabilized  ROMs  are 
observed  to  predict  solutions  that  are  in  excellent  agreement  with  the  CFD  results.  In  contrast  with  the  results  in 
Figs.  5  (bottom)  and  7,  adding  the  additional  amount  of  artificial  dissipation  is  effective  in  dissipating  the  unphysical 
instabilities  and  enables  the  ROMs  to  capture  the  original  CFD  solutions  very  well. 
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VI.  Conclusions 


Stability  issues  encountered  in  previous  studies  using  the  POD/Galerkin  method  for  reduced  order  model  (ROM) 
development  of  the  linearized  Euler  equations  are  investigated.  For  simplicity,  the  combustion  response  in  our 
model  equation  is  deactivated  and  only  the  gas  dynamics  are  calculated  to  obtain  the  POD  basis  vector,  which  are 
then  used  within  a  Galerkin  formulation  to  reduce  the  governing  partial  differential  equations  to  a  set  of  ordinary 
differential  equations.  Efforts  have  been  made  to  resolve  ROM  stability  issues  by  changing  and  scaling  the  solution 
variables  used  in  POD  mode  generation,  which  shows  some  success  but  still  does  notguarantee  a  universally  stable 
method.  Hence,  the  studies  in  this  paper  are  performed  on  two  major  aspects:  first,  to  understand  and  identify 
possible  causes  of  the  non-physical  instability  issues  present  in  the  ROM  development  for  vector  systems  of 
equations;  second,  to  test  and  evaluate  the  addition  of  extra  artificial  dissipation  to  mitigate  the  ROM  stability  issues. 

Even  though  only  simple  flow  dynamics  are  solved,  unstable  modes  are  still  encountered  in  the  resulting  ROMs 
and  the  inclusion  of  multiple  frequencies  in  the  forcing  function  appears  to  create  more  difficulties  in  obtaining  a 
stable  ROM.  Different  mean  flow  Mach  numbers  are  tested  to  generate  the  ROMs  and  it  is  found  that  the  ROM 
stability  issues  improve  as  the  Mach  number  gets  higher.  Moreover,  the  system  becomes  completely  stable  for 
supersonic  flows,  which  suggests  that  the  reduction  and  eventual  elimination  of  the  upstream  propagating 
characteristic  waves  is  beneficial  to  stability.  Following  that,  different  boundary  conditions  are  implemented  to 
control  the  amount  of  upstream  propagating  waves  present  in  the  CFD  solution.  It  has  been  shown  that  stable  ROMs 
can  be  generated  by  applying  either  non-reflecting  boundary  condition  downstream  or  Riemann  invariant  boundary 
conditions  both  upstream  and  downstream,  which  indicates  the  important  role  that  the  upstream  traveling 
characteristic  waves  have  on  the  overall  ROM  stability.  A  fundamental  insight  is  that  the  elimination  of  wave 
reflections  at  boundaries  essentially  reduces  the  Euler  equations  to  a  system  of  independent  characteristic  waves  that 
is  equivalent  to  multiple  scalar  wave  equations.  Indeed,  it  is  the  interaction  between  waves  that  arises  in  the  vector 
system  that  is  at  the  source  of  the  stability  issues. 

A  simple  solution  is  proposed  to  mitigate  the  instabilities  by  increasing  the  amount  of  artificial  dissipation  in  the 
the  Galerkin  projection  step.  The  additional  dissipation  is  shown  to  have  a  more  significant  effect  on  the  unstable 
modes,  while  they  do  not  affect  the  stable  modes  very  much.  Thus,  the  stabilized  ROMs  prove  to  be  able  to 
reproduce  the  original  CFD  solutions  with  excellent  agreement. 

Overall,  resolving  stability  issues  identified  in  the  ROM  development  for  vector  systems  of  equations  remains 
an  importantand  highly  challenging  task.  In  the  present  paper,  significant  progress  has  been  made  to  identify 
possible  causes  for  spurious  unstable  modes  in  the  ROM  construction  and  the  corresponding  stabilization  treatment. 
More  detailed  and  systematic  studies  are  needed  to  understand  the  underlying  causes  of  the  instability  and  to  insure 
that  the  proposed  addition  of  artificial  dissipation  remains  robust  for  more  complicated  problems. 
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